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Abstract .—The eastern gray squirrel, Sciurus carolinensis (EGS) has been introduced 
to California and has expanded its geographic range since initial introductions. In this 
study we projected the potential future geographic range of the EGS in California us¬ 
ing Maxent to create an ecological niche model. Location data were obtained over the 
time period of 2004-2015 from museum specimens, wildlife rehabilitation centers, the 
California Department of Public Health, the California Roadkill Observation System, 
and non-iNaturalist citizen science observations. Research grade data from iNatural- 
ist was obtained over the time period of 2004-2018. Range and habitat suitability 
maps were developed by mapping in ArcGIS. Three threshold selection methods were 
used to create different estimates of the potential future range of the EGS in Cali¬ 
fornia. The first method used the 10th percentile logistic threshold, the second used 
the minimum training presence logistic threshold, and the third used Jenks Natural 
Breaks. We propose that Jenks Natural Breaks has distinct advantages over the other 
two methods for estimating the potential future range of the introduced EGS in Cal¬ 
ifornia, because it provides information on the habitat suitability ranking throughout 
California, whereas the other methods only provide a binary suitable/unsuitable map. 


The objective of this study was to develop an ecological niche model (ENM) with an ap¬ 
propriate threshold value that could best identify potential range expansion of the invasive 
eastern gray squirrel, Sciurus carolinensis (EGS) within California and in the future, project 
potential areas of overlap with congeners. The EGS is native to the deciduous forests of 
the eastern United States (Koprowski 1994). The species was introduced to California in 
1939 at Stanford University and in 1943 at Golden Gate Park in San Francisco, but it 
may have been introduced earlier by settlers who frequently introduced species from their 
homes in the eastern United States (Byrne 1979). Populations of the EGS currently exist in 
developed and forested areas of California. They are currently widespread throughout cen¬ 
tral California, with concentrations around Sacramento, both peninsulas of San Francisco 
Bay, and San Jose, and smaller populations in the Central Valley. Additionally, populations 
are spreading from Santa Cruz into the Santa Cruz Mountains and Monterey Peninsula 
(Creley and Muchlinski 2017). 

Sciurus carolinensis carolinensis is most likely the major subspecies present in Califor¬ 
nia, as determined by the coat color and physical characteristics of observed squirrels. 
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S. c. carolinensis has a gray dorsum with a cinnamon wash sometimes present on the 
dorsum and hips. The tail is the same shade of gray, with a light white frosting on the 
tips of the hairs. A white eye ring is usually visible (Thorington et al. 2012). Some EGSs 
in California are melanistic, a common trait in the northern portion of the native range 
(Thorington et al. 2012). The EGS has a broad diet (Bertolino 2008), can establish a pop¬ 
ulation from a small number of founders (Wood et al. 2007), can survive and reproduce 
in urban, suburban, or natural habitats, and has a favorable public perception (Bertolino 
and Genovesi 2003), which all contribute to its invasive success. The species has been intro¬ 
duced to the western United States, Europe, Africa, and Australia (Bertolino 2008; Benson 
2013; Bertolino and Lurz 2013; Bertolino et al. 2014). Populations have been associated 
with negative ecosystem effects, the decline of native species, and damage to forests in the 
United Kingdom, Ireland, Italy, and parts of western North America (Gurnell et al. 2004; 
Bertolino 2008; Benson 2013; Bertolino and Lurz 2013). 

The original native range of the EGS consists of mature, continuous woodlands over 
40 ha in size, with diverse woody understories and tree species such as oak (Quercus), 
hickory (Carya), and walnut {Juglans) (Koprowski 1994). However, EGSs can also live 
in urban and suburban environments, even with relatively few mature trees (Thorington 
et al. 2012). EGSs move primarily along river corridors, secondarily on roads/right-of- 
ways, and thirdly on tracks/paths (Stevenson et al., 2013). Their ability to live in developed 
environments may significantly increase their dispersal capability. 

In order to understand the future potential distribution of the EGS, it is critical to model 
their potential landscape and ecological niche occupancy. Creley and Muchlinski (2017) 
mapped the species distribution within California as of 2015, but ENMs had not been 
made. The methodologies for creating ENMs for invasive species are not well established 
(Aguierre-Gutierrez et al. 2013; Uden et al. 2015). It was important to map a range of suit¬ 
able habitat estimates in order to prevent drawing conclusions from one arbitrary thresh¬ 
old, as cautioned against by Merow et al. (2013). Standard thresholds that have been used 
in previous studies to produce ENMs include the 10 th percentile logistic threshold (Belar- 
main Fandohan et al. 2015; Chalghaf et al. 2016; Sage et al. 2017), which could produce 
a very conservative estimate of potential future range for an invasive species, and the min¬ 
imum training presence logistic threshold (Beane et al. 2013; Coudrat and Nekaris 2013; 
Botero-Delgadillo et al. 2015), which could produce an overestimate of potential future 
range. We used these two standard thresholds in the present study, as well as a third method 
- selecting a threshold value based on Jenks Natural Breaks (Jenks 1967). Previously, Col- 
nar and Landis (2007) developed a regional risk assessment for the European green crab, 
Carcimis maenas , at Cherry Point, Washington, USA using Jenks Natural Breaks, and 
Schleier III and Sing (2008) used it to partition an overall risk score for the introduction of 
Gabusia affinis (western mosquitofish) into Montana watersheds. Beyond invasive species 
modeling, Jenks Natural Breaks have been used to classify groundwater into zones of vul¬ 
nerability for nitrogen contamination in Florida’s aquifers (Cui et al. 2016), to rank the 
susceptibility of locations to terrorist actions (Patterson and Apostolakis 2007), and to as¬ 
sess the risk of flooding in the Bengawan Solo River basin in Indonesia (Rahadianto et al. 
2015). 

The model used location and environmental data from the invaded range in California 
and the native range for S. carolinensis in the eastern United States, which encompasses 
a wider range of environmental conditions than those found within the current range in 
California, to estimate habitat suitability. Our methods may be applicable to further stud¬ 
ies on invasive species modeling, for which methods of estimating the potential range, as 
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opposed to the current range, are limited and have not been consistently evaluated 
(Aguierre-Gutierrez et al. 2013; Uden et al. 2015). 


Materials and Methods 

We obtained presence only location data for 2004 to 2018. The model included 3,725 
spatial location points of the EGS in California, as well as 8,988 points from across the 
native range in the eastern United States. We obtained presence points in California from 
iDigBio, the Global Biodiversity Information Facility (GBIF), Vertnet, wildlife rehabili¬ 
tation centers, the California Department of Public Health’s West Nile Virus Surveillance 
Program, the California Roadkill Observation System operated by the University of Cali¬ 
fornia, Davis (Waetjen and Shilling 2017), iNaturalist, and the authors. Location data from 
the native range of the EGS across the United States are from iNaturalist. We filtered the 
data to include only the native range of S. carolinensis , according to Koprowski (1994). We 
excluded regions from which the EGS is nonnative, and areas outside of the United States. 

VertNet, a National Science Foundation funded project, makes museum-curated 
biodiversity-data free and available on the web, while the Global Biodiversity Informa¬ 
tion Facility provides open international data. The iNaturalist sightings were filtered to 
include only those that were open access and research grade, which include an observation 
date, photo, coordinates, and in which the species identification has been verified by at least 
one other user. Records in biodiversity databases are constantly changed and updated, so 
all data from iDigBio, GBIF, and Vertnet were obtained on 31 August 2015. Records from 
iNaturalist were obtained on 12 October 2018. Some redundancy may exist between the 
databases, but Maxent automatically removes replicates, so these duplicate observations 
did not change the projections. 

Since all sources are likely to include some misidentifications of related species identified 
as EGS, reports from outside of the previously published range were scrutinized for accu¬ 
racy. Field surveys were conducted in regions that had not been included in prior range 
maps but that had numerous reports, including the Santa Cruz Mountains, the Central 
Valley, and southern California. Records that could not be corroborated were expunged. 
In the case of the California Roadkill Observation System data, 82% of EGS with images 
were correctly identified. The remainders were misidentified as California ground squirrel, 
Douglas squirrel, or western gray squirrel (WGS). 

Data used in this manuscript are available for use by others under a Creative Com¬ 
mons By Attribution Non-Commercial 4.0 International License. Observations ob¬ 
tained through iNaturalist are utilized under a Creative Commons By Attribution Non- 
Commercial License or from observations that are in the Public Domain. The names of 
GBIF and iNaturalist contributors can be obtained through a search of the California 
location data using the data posted at doi: 10.13140/RG.2.2.24275.84004 and the native 
range of the United States data at doi: 10.13140/RG.2.2.17564.95360. 

We converted location data into geographic coordinates with Google Maps. We spa¬ 
tially rarefied the data to one point per 51.8 ha (0.25 mi 2 ) using the Spatially Rarefy tool 
in the SDM toolbox (Brown 2014) in order to eliminate spatial clusters that could cause 
the model to be overfit to the environmental biases of those points (Boria et al. 2014). The 
rarefied presence data were reduced from 12,713 to 5,627 points. Since the predictive accu¬ 
racy of the maps and the ability of the models to project presence data was more important 
than identifying the tolerance ranges of the species, we did not remove highly correlated 
variables (Merow et al. 2013). 
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For the environmental background we used the bio 1 through bio 19 variables from BIO- 
CLIM1 remote sensing data (Hijmans et al. 2005), which include quarterly and annual 
temperature and precipitation trends, with other annual environmental trends. Addition¬ 
ally, we used monthly precipitation, monthly maximum temperature, monthly minimum 
temperature, altitude, impervious surface, land cover, and tree canopy from the United 
States Geological Survey (USGS 2016). The full list of biotic and abiotic factors used in 
our model is available at doi: 10.13140/RG.2.2.15802.70085. We paired the location pres¬ 
ence data with environmental background data throughout the contiguous United States. 
The model can predict habitat suitability for areas that the species is equally likely to reach 
(Merow et al. 2013). The environmental background data covers all reasonable possible 
distributions (Saupe et al. 2012). 

We selected Maxent to create the ENM because of its high performance at estimating 
local occurrences with small to medium sample sizes (Elith et al. 2006; Aguierre-Gutierrez 
et al. 2013; Ng and Jorda, 2001). In order to allow the model to reach the default 0.00001 
convergence level, we allowed a maximum of 5,000 iterations (Young et al. 2011), set the 
regularization parameter to the default value of one in order to reduce overfitting (Merow 
et al. 2013), and created 15 replicates for the model using subsampling (Young et al. 2011). 
We set aside 25% of the data for testing, and used 75% for training the model. We used the 
random seed option in order to increase the randomness of the runs (Jobe and Zank 2006; 
Young et al. 2011). We adjusted the sample radius to -100, and did not extrapolate. The 
logistic output using the default t value of 0.5 was selected because the actual probability 
of the EGS being present in suitable habitat is unknown. 

We produced habitat suitability maps by importing the rasters for the average of the 
fifteen replicates for each model into ESRFs ArcMap 10.3.1 (ArcGIS® 10.3.1; Esri soft¬ 
ware) using the NAD 1983 California (Teale) Albers (Meters) projected coordinate system 
and the GCS North American 1983 (NAD 1983) datum. When coordinates were provided 
without a datum, they were assumed to be in NAD83. Each pixel encompasses approxi¬ 
mately 720 square meters. We clipped the output rasters to the shape of California, using 
the United States Census Bureau’s Tiger/Line 2010 (United States Census Bureau 2010). 
We used topographic and political basemap layers from ESRI. 

We used three thresholds to create potential range estimates. The first estimate of the po¬ 
tential range was established by using the 10th percentile training presence logistic thresh¬ 
old, which set the threshold at the level where 90 percent of the presence points were in 
raster squares with at least the threshold score. We created a second estimate of the po¬ 
tential range using the minimum training presence logistic threshold, which is set at the 
level of the lowest scoring occupied square. We manually classified the resulting raster of 
California for the first two thresholds by setting the upper bound of the unsuitable habi¬ 
tat class at the threshold level for each map, as indicated in the Results table created by 
Maxent. 

Third, we used Jenks Natural Breaks to group the resulting raster into five classes 
based on natural breaks in the data. The approach grouped the relative habitat suitabil¬ 
ity ranking of each raster square (Rahadianto et al. 2015) by similar values and max¬ 
imized the difference between the classes (ESRI 2016). We created 5 habitat categories 
(HC) from the groupings of data. HC1 represents where the species is now found and 
adjacent highly suitable habitat, HC2, HC3, and HC4 represent decreasingly suitable, 
but still suitable habitat, while HC5 represents unsuitable habitat. Therefore, the thresh¬ 
old for suitable habitat was set at the lowest relative habitat suitability ranking score 
of HC4. 
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Results 

The high value of the Test Area Under the Curve (AUC) of the Receiver Operating Char¬ 
acteristic (ROC) Plot (0.8272 on a scale of 0.5 = random association of data to model, to 
1.0 = perfect association of data to the model) and Training AUC (0.828) indicate an 
excellent fit of the data to the model produced by Maxent. The plot is available at doi: 
10.13140/RG.2.2.12531.78883. However, since the ROC Plot is constructed from all pos¬ 
sible threshold values and the plot does not yield information about any specific threshold 
value, the selection of a biologically meaningful threshold value for graphing results of the 
model is critical for obtaining the most accurate map of potential future range. 

Three different threshold values were used for the range of suitable to unsuitable habitat 
found using the ENM for the EGS in California. Use of the 10 th percentile logistic thresh¬ 
old (Fig. 1) produced the most limited estimate of potential future range. In the binary 
map, suitable habitat scored above the threshold value of 0.3051 (Table 1). The potential 
range produced using this threshold value is very similar to the current range of the species 
in northern and central California as of 2018. Therefore, the time frame of projected range 
expansion is very limited. EGSs currently exist outside of habitat that is predicted suitable 
through use of this threshold, as the method excludes the 10 percent of observation lo¬ 
cations with the lowest relative scores in order to produce the threshold value. EGSs are 
currently found outside of projected suitable habitat near Salinas, Tracy, Modesto, and east 
and south of Sacramento along the Sierra Nevada foothills from north of the American 
River to Columbia. 

Use of the minimum training presence logistic threshold (Fig. 2) also produced a binary 
output map, but in this case the logistic threshold value was reduced to 0.0299 (Table 1). A 
very large increase in projected suitable habitat for the EGS in California results from the 
reduction in threshold value. The potential range using this threshold value includes most 
of the state, with the exception of the Mojave Desert, the very north-central and northeast 
portion of the state, the higher elevations of the Sierra Nevada Mountain Range, and a 
southern portion of the San Joaquin Valley. 

The Jenks Natural Breaks method produced a binary map with suitable habitat above 
the threshold of 0.0530, but the specificity of information within the suitable habitat cat¬ 
egory was much greater than in either of the previous two methods (Tables 1, 2). In the 
map classified using Jenks Natural Breaks (Fig. 3), the two highest relative habitat cate¬ 
gory rankings (HCl and HC2) are centered on the San Francisco Bay area, south through 
the Santa Cruz Mountains to habitats on the Monterey Bay Peninsula, east through 
Sacramento into the foothills of the Sierra Nevada Mountains, and along the southern 
California coasts. Concentric areas of decreasingly suitable habitat (HC3 and HC4) sur¬ 
round the most suitable areas and continue south along most of coastal California and the 
coastal/inland mountain ranges to the border with Mexico. HC3 and HC4 are also found 
along river corridors exiting the Sierra Nevada Mountains. Unsuitable habitats (HC5) are 
in and east of the Sierra Nevada Mountains, the northern Cascade Range, a major portion 
of the San Joaquin Valley south of Sacramento, and the deserts of southern California. 

Discussion 

Our results indicate that the method for selecting a threshold value to establish a break 
between suitable and unsuitable habitat with Maxent is an especially critical issue for an in¬ 
troduced species which is expanding its geographic range. We have shown that a high AUC 
value by itself is not sufficient to support the predictive accuracy of a single map based 
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Fig. 1. Estimate of suitable habitat predicted by the model using the 10th percentile training presence 
logistic threshold based upon eastern gray squirrel locations from 2004 to 2018. 

upon a single arbitrarily selected threshold value. AUC is not a perfect, objective measure 
of the predictive power of the model, but few alternatives are available for presence only 
data (Merow et al. 2013). The commonly used 10th percentile training presence logistic 
threshold as well as the minimum training presence logistic threshold are clearly arbitrary 
values with little biological basis for selection. And, as shown by the maps presented in this 
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Fig. 2. Estimate of suitable habitat predicted by the model using the minimum training presence logistic 
threshold based upon eastern gray squirrel locations from 2004 to 2018. 


paper, these thresholds produce vastly different projections of suitable habitat. The use of 
Jenks Breaks to establish a threshold value is based upon the natural clustering of values 
in the logistic output of the model and hence should be less arbitrary, and more meaning¬ 
ful, in terms of the biology of the species being studied. This method also provides graded 
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Table 1. Threshold breakpoints. 



10th percentile training 

Minimum training presence 

Jenks 

Class 

presence logistic threshold 

logistic threshold 

Natural Breaks 

Suitable 

0.3051 - 1 

.0299-1 

0.0530 - 1 

Unsuitable 

0- 0.3051 

0-0.0299 

0 - 0.0530 


“likelihoods” of predicted occupancy, which may be more easily tested in the future with 
new occupancy records. 

The 10th percentile training presence logistic threshold provided the most conservative 
estimate of the potential range by assuming that some of the presence data may be misiden- 
tified, improperly reported, or outside of the area in which the EGS can persist, and then 
removing those points (Young et al. 2011; Uden et al. 2015). The map is inherently skewed 
toward the realized niche of the EGS in California, as the species has only been introduced 
to a few locations within the state. It is closely aligned with the original data set and may 
be biased toward human accessible areas. The threshold may be too conservative, as the 
EGS is still spreading, and is likely to tolerate habitat with conditions at least as extreme 
as those in the current range in California. The EGS already inhabits areas outside of the 
region projected as suitable in this map. 

The map created with the minimum training presence logistic threshold reflects a much 
broader range of conditions throughout California and is our least conservative estimate of 
the future range. However, use of this threshold value most likely over-predicts the poten¬ 
tial range because the inclusion of a single erroneous location point within the data could 
greatly affect the map. The map using this lower threshold provides an accurate estimate 
of potential range only if it is absolutely certain that all location data have been correctly 
identified as EGSs, and all individuals and populations currently exist in suitable habitat. 
While it is highly likely that the EGS will expand its range to include some of the areas 
mapped by this method, it is highly unlikely that the species will inhabit all of the area. If 
the EGS is introduced to or expands it range to new areas and survives, the distribution 
is likely to more closely resemble the maximum estimate map than the minimum estimate 
using the 10th percentile logistic threshold. 

The map created with Jenks Natural Breaks provides the most useful, detailed classifica¬ 
tion of the relative habitat suitability rankings, and therefore the risk of invasion to various 
parts of California. We assume that the lowest ranking habitat, HC5, is unsuitable based 
upon the lowest grouping of values (0.000 to 0.0530), whereas each higher rank represents 
increasingly suitable habitat. As the EGS continues to spread, the habitat suitability rank¬ 
ing may increase in areas with habitat similar to newly invaded areas and therefore we 
cannot say that this map is a permanent ranking of habitat suitability. For example, areas 


Table 2. Habitat classes created with Jenks Natural Breaks. 


Habitat class 


Relative habitat suitability ranking 


1 0.4627-0.7515 

2 0.2917 -0.4627 


3 

4 

5 


0.1562-0.2917 
0.0530 -0.1562 
0 - 0.0530 
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Fig. 3. Estimate of suitable habitat predited by the model using Jenks Natural Breaks, with five habitat 
classes (HC), based upon eastern gray squirrel locations from 2004 to 2018. HC1 represents the most suit¬ 
able habitat. HC2-HC4 represent decreasingly suitable habitat, and HC5 is considered unsuitable habitat. 

now ranked as HC3 could in the future be ranked as HC2 if many populations of EGSs 
become established within HC3. We believe that setting the threshold for suitable habitat 
above HC5 created the most reasonable estimate of the potential future range by including 
a broader range of habitat than the minimum estimate, but excluding the most extreme 
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desert, mountain, and Central Valley habitats of the maximum estimate. The Jenks Nat¬ 
ural Breaks map is, at this time, our most reasonable estimate of the future range of the 
EGS in California. 

We selected Maxent to create the ENM because of its high performance at estimating 
local occurrences. Absence data were not available or informative, so a presence-absence 
experiment could not be done (Yackulic et al. 2012). Additionally, presence only data may 
be better than presence-absence data for estimating fundamental niches because invasive 
species have not yet inhabited their full potential range (Jimenez-Valverde et al. 2008). For 
presence only data, machine learning methods such as Maxent consistently outperformed 
earlier methods, such as Bioclim or regression models in predictive success (Elith et al. 
2006). It is especially good for small and medium data sets because the generative learning 
method uses an algorithm to build a model, as opposed to a discriminative model that 
estimates the value of one categorical variable based on the other (Ng and Jordan 2001; 
Aguierre-Gutierrez et al. 2013). 

Location data were limited to 2004-2018 because habitat loss has been extensive in Cal¬ 
ifornia over the past century, to ensure the use of current species distributions, and be¬ 
cause satellite collection for BIOCLIM only started in 1972 (United States Geologic Sur¬ 
vey [USGS] 2015). Overall, 2,825 of the 3,725 EGS records in California, 75.84%, were 
obtained from wildlife rehabilitation records. Records indicated the location that the per¬ 
son submitting the squirrel reported finding it. Rehabilitation facilities are very likely to 
accurately identify the species, but may cause geographic bias toward the residential areas 
surrounding each center. Citizen science data offers the benefit of having broad sources, 
but may be skewed toward urbanized areas, roads, or other easily studied sites (Baldwin 
2009). The data are not random, but rarefaction has been used to equalize the sampling 
effort among areas, which allowed us to make a reasonable inference of the species distri¬ 
bution (Yackulic et al. 2012). The assumption that detection probability is constant across 
sites need not be met because citizen science and historic records simulate repeat-visits to 
each site (Yackulic et al. 2012). 

The higher elevations of the Sierra Nevada were classified as unsuitable using all thresh¬ 
olds, but elevation may have been the only factor causing the unsuitable classification, as 
the highest point in the native range of S. carolinensis is at 6,643 feet, at Clingman’s Dome 
in Tennessee, which is the highest point in the state and the Smoky Mountains. Elevation 
alone may not actually exclude the EGS. The species’ invasive ability worldwide, especially 
into areas with Mediterranean climates, such as South Africa and Italy (Gurnell et al. 2004; 
Bertolino 2008; Benson 2013; Bertolino and Lurz 2013) suggests that it could acclimate to 
conditions in much of California. 

Finally, with future studies the potential invasion by the EGS into occupied and un¬ 
occupied eastern fox squirrel (EFS) or WGS habitat is an excellent test case for several 
hypotheses regarding invasion ecology. These hypotheses include, 1) “biotic resistance”, 
which suggests that high-biodiversity ecosystems are more resistant to invasion than low- 
diversity ecosystems; 2) “enemy release”, which posits that the absence of enemies (i.e., 
competitors and predators) increases the likelihood of invasion; and 3) “propagule pres¬ 
sure/introduction effort”, which proposes that the introduced population size and the fre¬ 
quency of introduction can contribute to successful invasion ( Jeschke 2014). The first hy¬ 
pothesis could be tested with the EFS by comparing the rate and success of invasion from 
neighboring areas of low and high-biodiversity systems (e.g., oak woodlands in the Sierra 
Nevada mountain foothills). The second hypothesis could be tested by comparing EGS 
invasion and persistence success with the presence of competitors (e.g., EFS or WGS) 
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and predators, such as hawks and owls. Because the EGS tolerates urban and suburban 
areas, predators may be less prevalent in areas where they have succeeded. The third hy¬ 
pothesis could be tested using a combination of estimates (e.g., from historical accounts) 
or measurements (e.g., for newly-discovered populations) of founder population size and 
continued introduction (e.g., from connected populations) and reproductive rate to predict 
expansion and persistence. There is already some evidence that the EGS may be replac¬ 
ing the EFS on the western side of the San Francisco Bay (Creley and Muchlinski 2017) 
so assessing potential causal factors for replacement is important. Finally, the model uses 
current environmental conditions, and could be adapted to include estimates for changes 
in climate, which may reduce the amount of suitable habit, shift it to higher elevations, or 
render hotter and more extreme areas of currently suitable habitat unsuitable in the future. 
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